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Abstract. Event horizons are the denning physical features of black hole 
spacetimes, and are of considerable interest in studying black hole dynamics. 
Here, we reconsider three techniques to find event horizons in numerical 
spacetimes: integrating geodesies, integrating a surface, and integrating a level- 
set of surfaces over a volume. We implement the first two techniques and find 
that straightforward integration of geodesies backward in time is most robust. 
We find that the exponential rate of approach of a null surface towards the event 
horizon of a spinning black hole equals the surface gravity of the black hole. In 
head-on mergers we are able to track quasi-normal ringing of the merged black 
hole through seven oscillations, covering a dynamic range of about 10 . Both at 
late times (when the final black hole has settled down) and at early times (before 
the merger), the apparent horizon is found to be an excellent approximation of 
the event horizon. In the head-on binary black hole merger, only some of the 
future null generators of the horizon are found to start from past null infinity; 
the others approach the event horizons of the individual black holes at times far 
before merger. 



PACS numbers: 04.25. Dm, 04.30.Db, 04.70.Bw 

1. Introduction 

The two body problem in general relativity has been the focus of extensive work for 
many years and, because there is no analytic solution, it must be solved numerically. 
Binary black hole mergers are expected to be one of the most astrophysicaUy common 
sources of gravitational radiation for detectors such as LIGO p] [2] • Recent advances 
in simulating binary black hole mergers include the development of the generalized 
harmonic evolution system [3] and the moving punctures technique [H [5]. In the 
last several years the field has reached a stage where binary black hole simulations 
are becoming routine. Numerical simulations have been remarkably successful in 
expanding our understanding of binary black holes, but challenges remain. 

One particular challenge is to be able to more accurately locate the holes during 
the merger. There are two useful concepts to describe the location of black holes in a 
spacetime, apparent horizons (AH) and event horizons (EH). An EH is the true surface 
of a black hole: it is defined as the boundary of the region of the spacetime that is 
causally connected to future null infinity. Because the definition of the EH involves 
global properties of the spacetime, one must know the full future evolution of the 
spacetime before the EH can be determined exactly. This difficulty has led researchers 
to instead identify black holes with apparent horizons, which are defined in terms of 
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the expansion of null congruences, [g Indeed, AH finders are highly developed and 
have been the subject of extensive work (see, e.g. the review [5]) Unlike an EH, an 
AH can be located from data on a single spacelike hypersurface, i.e. on each timestep 
of a numerical evolution, without knowing the future evolution of the spacetime. The 
AH is often an effective substitute for the EH for several reasons. First, according to 
the cosmic censorship conjecture, if an AH is present, it must be surrounded by an 
EH. Second, if an AH is present on a spacelike hypersurface through a stationary 
spacetime, it coincides with the EH. Finally, in numerical simulations, apparent 
horizons generally show behaviour attributed to event horizons: For instance, the 
area of the AH typically does not decrease and it is usually almost constant whenever 
the spacetime is only mildly dynamic. In fact, apparent horizons have motivated the 
development of "isolated" and "dynamical" horizons (see [7] for a review). These 
surfaces satisfy analogues of the laws of black hole thermodynamics, although they 
are defined quasi-locally, rather than globally. 

However, using the AH to locate the holes is not always appropriate. For instance, 
the AH is slicing dependent, while the EH is not. Indeed, the Schwarzschild spacetime 
can be sliced in such a way that no AH exists [8]. Furthermore, even on slicings on 
which an AH is present, there are few precise mathematical statements about how 
"close" AH and EH are. Finally, AH and EH behave qualitatively differently during 
a black hole merger: The EH |§| around each black hole expands continuously until 
the two components of the EH join into one, whereas a common apparent horizon 
appears discontinuously quite some time after the EHs have merged. The common 
AH encompasses the two individual AHs, which continue to exist as surfaces of zero 
outgoing null expansion for some time after the merger. 

Early EH finders [HI QH] followed null geodesies forward in time and determined 
whether or not each geodesic eventually escapes to infinity. Following geodesies 
forward in time is unstable in that slightly perturbed geodesies will diverge from 
the EH and either escape to infinity or fall into the singularity. Furthermore, a large 
number of geodesies with different directions must be sampled at each point and at 
each time step to determine if one of these succeeds in escaping to infinity [TO]. To 
reduce the number of sampling points, the EH search in [lOj was performed on a series 
of time slices proceeding backward from late to early times; to find the EH on each 
time slice, they integrated geodesies forward in time, using the already-located EH at 
the later time as an initial guess. 

Since outgoing null geodesies diverge from the event horizon when going forward 
in time, when going backward in time they will converge onto the event horizon |llU12j . 
All recent EH finders use this observation, and follow null geodesies or null surfaces 
backward in time El H2 H3 HI 113 HSl H3 HH1 US] ■ 

Several algorithms have been developed to follow null geodesies backward in time. 
These can be divided into three types, which we shall refer to as the "geodesic method" , 
the "surface method" and the "level-set method" . The geodesic method works by 
simply integrating the geodesic equation, as done by Libson et. al. [P2j . Libson et. al. 
express concerns that the geodesic method may be susceptible to tangential "drifting" 
of the geodesies. However, this is not evident when the method is applied to the science 
applications in that paper, nor do we find tangential drifting in our simulations. To 
avoid any issues with drifting, Libson et al. introduced the surface method: a complete 

X More precisely, we define AH as the outermost marginally outer-trapped surface, where an outer- 
trapped surface is a topological 2-sphere with zero expansion along outgoing null normals. 
§ More precisely, the 2-surface formed by the intersection of the spatial slice and the EH 
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null surface (rather than individual geodesies) is evolved backward in time. In [TTJ [T2] 
this surface was parameterized based on axisymmetry (although the parameterization 
of [111 112] cannot handle generic axisymmetric situation, cf. Section |2 . 21 below) . and 
many interesting results on the structure of caustics and the geometry of the horizon 
for axisymmetric spacetimes were obtained in [13l [15]. Diener [19] and Caveny et. 
al [HI [T71 [18] independently introduced the level-set method by recasting the surface 
method in a way that does not assume symmetry: rather than evolving a single 2- 
D surface, they evolve a volume-filling series of surfaces given as the level-sets of a 
spacetime function f(t, x l ). To avoid exponentially steepening gradients of /, Caveny 
et al introduce an artificial diffusive term, whereas Diener reinitializes / whenever 
necessary. 

This paper re-examines these techniques for event horizon finding in the context of 
the Caltech/Cornell Spectral Einstein Code (SpEC), which provides an infrastructure 
for highly accurate simulations of Einstein's equations for single and binary black holes. 
Recent work includes highly accurate computations of gravitational waveforms from 
inspiraling binaries [20] [21j [22] . The availability of high accuracy binary evolutions 
motivates the development of very precise event horizon finding techniques in order to 
extract all possible physics from these simulations. Therefore, this paper reconsiders 
the three techniques mentioned above in the context of general binary black hole 
mergers without any symmetries. 

We implement the geodesic method, and generalize the surface method to 
arbitrary situations without symmetries. Both methods are then applied to single 
Kerr black holes, and a head-on binary black hole merger. In both cases, the geodesic 
method is found to be more robust. We encounter two fundamental problems with 
the level-set method, and therefore halted our efforts to implement it in SpEC. 

This paper is organized as follows: In Section [2] we explain the three methods 
in more detail and give details of our numerical implementation. Section [3] presents 
results for a single Kerr black hole, and in Section [4] we apply the techniques to a 
head-on BBH merger, where we extract ringdown behaviour and the behaviour of the 
individual event horizons before merger. We close with a conclusion in Section [5] 

2. Methods 

All EH-finding techniques considered here proceed backward in time and must 
therefore be performed after the numerical evolution of the spacetime has been 
completed. We assume that we have access to the spacetime metric in a 3+1 
decomposition 

ds 2 = -N 2 dt 2 + jij (dx l + P l dt){dx 3 + P>db), (1) 

where N is the lapse, /3* is the shift, and %j is the 3-metric on the slice. Latin 
indices = 1,2,3 denote spatial dimensions; below we will use Greek indices 

to denote spacetime dimensions, a, j3 . . . = 0, 1, 2, 3. The time t in (JTJ) represents the 
coordinate time of the numerical evolution. Typically, the metric data jij, (J 1 , and 
N are available at discrete times and at discrete spatial grid points. Evaluating the 
values of the metric components elsewhere requires interpolation. 

A black-hole merger exhibits several characteristic features of relevance to EH 
finders, as illustrated in Figure □U(cf. [H HI QH). At times sufficiently far prior 

|| While Figure [T] is meant as an illustration, it presents actual data from the head-on binary black 
hole merger discussed in Section [4] The time given at each frame of Figure [T] will aid in the discussion 
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Figure 1. Cross-sections through event and apparent horizons during a BBH 
merger. Before the merger, t < tcEH > the surface includes the set of generators 
that will merge onto the event horizons through the cusps in the individual event 
horizons (green dashed curves). The point C is the point of symmetry for the 
head-on merger, which will be used in Section 12.31 




to merger, the EH and AH are expected to coincide closely (and indeed, we confirm 
this below for our simulation). The green dashed curves in Figure [T] represent future 
generators of the event horizon, i.e. null geodesies that will merge onto the event 
horizon through cusps in the individual event horizons. These cusps are clearly visible 
at time t — 13. 5M where the individual EHs have diverged significantly from their 
respective AHs. At <ceh = 14.6M the two previously disjoint components of the event 
horizon join. We shall refer to this time tcEH as the merger of the black hole binary. 
After the merger, the event horizon of the merged black hole can be seen relaxing 
towards its final time-independent shape. The common apparent horizon appears 
at icAH = 17.8M, and approaches the event horizon as the evolution proceeds; at 
t = 80M, the AH coincides almost exactly with the event horizon. 



2.1. Geodesic Method 

The most straightforward way to follow light rays is to simply integrate the geodesic 
equation [SEH QH E], 

dA 2 a ^~dX dX ~ ( ' 
where = (^(A) is the position of the photon on the geodesic, parameterized by an 
affinc parameter A, and T^o are the spacetime Christoffel symbols. 

Since we have access to our spacetime as a function of the evolution time 
coordinate t, it is convenient to rewrite ([2]), replacing A by t along the geodesic. 



in Section [4] 



Revisiting Event Horizon Finders 



5 



Writing q^ = dq^/dt, and a — dX/dt, we find: 
dq" 1 



m =-^» (3) 

(4) 



dA 2 

Substituting into the geodesic equation we get 



r = -r-K q a q?. (5) 



a 

,0 



The quantity a is determined by the requirement that q = t, i.e. that at parameter 
value t along the geodesic, the geodesic is on the corresponding t — const hypersurface 
of the evolution. This implies qf* = [l,<f] and = [0, q 1 ]. Setting p = in ([5]) gives 
2 = r^g"^. The spatial components of ([5]) are the desired evolution equation for 
the spatial coordinates as a function of coordinate time, 

q l ^T%q a q"q l -K q a q - (6) 

We convert this set of ordinary differential equations to first order form by defining 
p l = q l , which gives 

g*=P*, (7«) 

f = r>VV - r>V, (76) 

facilitates the use of standard ODE integrators like Runge-Kutta methods [UJ 121] • 

While integrating geodesies is not new [TUHH] , re-expressing the geodesic equation 
in terms of coordinate time seems to be new. It appears that the primary reason this 
technique has been phased out in favour of the two techniques described below is the 
concern that, in a full 3D implementation, slight tangential velocities may be imparted 
to the outgoing null geodesies through numerical inaccuracies, and that this tangential 
drift of geodesies could result in unphysical caustics. These concerns are discussed in 
detail in [T2|, where the idea of representing the whole surface, rather than individual 
geodesies, was introduced. This was justified on the basis that for a surface, tangential 
drift is irrelevant. However, while it is possible that tangential drift can be significant 
for very coarse, low-resolution simulations, we see no evidence that tangential drift 
affects our numerical tests of the geodesic method. 

We finally like to point out that if one evolves pi — g^p^ instead of p l (cf. 
[TO]), then the evolution equations depend only on spatial derivatives of the spacetime 
metric. Evolving pi therefore results in computational savings, because the time 
derivatives of the metric need not be stored or interpolated. This will be investigated 
in a future work. 



2.2. Surface Method 

The idea of the surface method dates back to Libson et al. [T2], who used it in 
axisymmetry. The goal is to evolve a 2-dimensional surface St backward in time 
such that it traces out a null hypersurface Af. The time coordinate t is inherited from 
the black hole simulation for which event horizons are to be determined, i.e. St is the 
intersection of Af with the spatial hypersurfaces S t of the evolution, as indicated in 
Figured Before the black hole merger, t < tcEH, the surface St consists not only of 
the two disjoint parts of the event horizon, but also includes the future generators, 
which are indicated by the green dashed curves in Figure [TJ The union of these three 
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Figure 2. A slice of the event horizon (St), produced by the intersection of a 
spatial hypersurface £t with the world tube of the event horizon J\f , showing n M 
(the timelike normal to £(), s M (the spatial normal of St), and £ M (the null normal 
to the EH world-tube AO- 



components is a smooth self-intersecting surface with the topology of a sphere (as 
suggested by Kip Thorne fTTl IT2"] ). 

Let us first consider how to represent the surface to be evolved. Apparent horizon 
finders often parameterize a surface by giving the radius, relative to a fixed point, as 
a function of angular coordinates, i.e. r = f(9,<p). Such a star-shaped surface is 
insufficient here, because the surface will be self-intersecting for t < tcEH and will 
cease to be star-shaped even before then (see Fig. 1.1 of [6]). The axisymmetric EH 
finder presented in [12] parameterized the surface by p — s(z, t), where z is a coordinate 
along the axis of symmetry, and p is the cylindrical radius. This allows some mild form 
of self-intersection, like, for instance, the t — 13. 5M snapshot in Figure [TJ However, 
at earlier times, the locus of future null generators of the horizon "bulges outward" 
and becomes multivalued when considered as a function of z, cf. t = 9M in Figure [T] 
In this case, the parameterization of [12] fails even for an axisymmetric configuration. 
In this paper, we use a parametric representation of St, i.e. r l — r l (t, u,v)\ . The full 
3-dimensional null hypersurface Af being constructed is represented as a 3-parameter 
surface in spacetime: 

r»(t,u,v)= [t,r%u,v)]. (8) 

We wish to find an equation that will allow us to evolve St in such a way as 
to trace out the null 3-surface Af '■ Further, we would like this equation to have the 
property that for fixed (uq, vq), the curve r^(t, uq, Vq) traces out a null geodesic. This 
allows us to directly compare the surface obtained by the surface method to the surface 
obtained for equivalent initial conditions by the geodesic method. 

For the curve ^{t, u,v)\ uv to be null, its tangent <9r M (t, u, v)/dt must be outgoing 
and null, i.e. 

dr^ 

(9) 

where is a null normal to St- can be written as 

£1" = C {7l" + S''), (10) 

where n M is the timelike unit normal to T, tl s M is the spatial outward-pointing unit 
normal to St (cf. Figure [2]) and c is an overall scaling. Consistency of J8} and ([9]) 
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requires that is normalized such that £* = 1. To find the value of c from the 
condition l l — 1, first notice that from the 3+1 decomposition, 

»" = -^[l, -/?*], (11) 

where N and /3 1 are the lapse and shift fields. Also, since s M lies within the spatial 
slice Sj, we may write s M = [0, s 1 ] , so that (JTDJ) becomes 



2V' iV^ 



(12) 



Thus £ — 1 implies c = N, and we can write our final evolution equation for the 
spatial components of r l , 

dr % 

— =Ns*-(3\ (13) 

In order to find the unit normal s 1 to the spatial surface S t , we follow the standard 
procedure for a surface parameterized as r l (u, v), i.e. 

■7 dr-i dr^ 

r^'toW (14a) 

P = y/lij^, (146) 

s^p" 1 ^. (14c) 

where ey^ is the antisymmetric tensor and where we have chosen the sign of the root 
such that s 1 points outward for a right-handed choice of coordinates. 

This evolution equation for the surface method ([Us]) is very different from the 
evolution equations for the geodesic equation (f7b|) - ([7a|) . The surface method does 
not require derivatives of the metric, but derivatives d u r l , d v r l along the surface; 
the geodesic method, in contrast, requires derivatives of the metric, but treats each 
geodesic completely independently. Nevertheless, due to our choice of evolution 
equation @, each point on the parameterization of the surface traces its own geodesic; 
see |Appcndix A| for a proof. 



2.3. Level-set Method 

The level-set method QH HI 021 El [HI OH QJ] utilizes a function / = /(t,^) defined 
on the full spacetime (or at least, a region of spacetime covering the vicinity of the 
expected location of the EH). The function / is determined such that / = const 
contours (i.e. level-sets) represent null surfaces i.e. g al3 d a fd/3f — 0. In the 3+1 
decomposition, this becomes [T71 [TSJ [TH] , 

<h.f r»,f : A\ "'-''V, /(',/. (15) 

where the + accommodates both ingoing and outgoing null surfaces, with the minus 
sign being appropriate for outgoing null surfaces if the gradient dif is outward- 
pointing. 

Libson et al [H] had previously made use of (fTS"]) . but parameterized the / = 
contour based on axisymmetry. The motivation of evolving (|15[) directly in the volume 
is to remove any assumptions of symmetries. 

Unfortunately, when trying to implement the level-set method in SpEC, we 
encountered two fundamental problems. The first difficulty is related to the 
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characteristic speed of the level-set method. Simply put, all / = const contours 
approach the event horizon, therefore new contours need to be filled in at the 
boundaries of the region in which / is evolved (i.e. the outer boundary and possibly 
one or more inner boundaries if black hole excision is employed). To see this, note 
that the characteristic speed of (|T5)) relative to a spatial direction Hj is 



i, = Nf,: — d f (16) 

where the sign of the first term depends on the gradient dif being outward pointing. 
For most coordinate systems of interest, lapse N and shift [3 l behave such that v > at 
the outer boundary and at any excision boundaries (if present). When integrating (|15[) 
backward in time, well-posedness requires boundary conditions at these boundaries. 
Our preferred numerical techniques are spectral methods because of their promise 
to achieve exponential convergence for smooth problems. Spectral methods are 
very sensitive to the existence of an underlying well-posed continuum problem and 
therefore require boundary conditions. Unfortunately there is no particular physical 
reasoning to suggest a choice of boundary condition. While essentially any choice of 
boundary condition that results in / being continuous rendered our spectral level-set 
implementation stable, and convergent to at least first order, we have been unable 
to find a boundary condition that ensures that / remains smooth and thus leads to 
the desired exponential convergence, not even in the single black hole case. A full 
finite-difference evolution of / would be less sensitive to the lack of proper boundary 
conditions (see [12]), but would be much slower for finding an EH in spectral-code 
metric data (due to interpolations from the spectral to the finite-difference grid) and 
much less accurate. 

The second fundamental difficulty lies in singular behaviour of the function / in 
certain cases. Let us consider an equal- mass head-on merger as depicted in Figure [T] 
Assume / to be smooth, and let us focus on the value of / at the point of symmetry, 
marked with C in Figure [TJ We assume that dif is outward-pointing near the event 
horizon. At late times, after the merger, / will be negative at C, because C is inside 
the event horizon. Throughout the whole simulation, dif = at C by symmetry, and 
therefore, (|15p implies that dtf = there, so that / at C remains fixed at a finite 
negative value. At merger, however, the / = contour passes through C. Therefore, 
/ must be singula^. Any method for solving the level-set equations that assumes a 
smooth and regular solution (including finite-difference methods that do not explicitly 
treat the singularity) will therefore produce results that differ from the exact solution 
at the singular point. In [19], one-sided finite-difference stencils are carefully chosen 
so as to not differentiate across the singularity. 

Because of these two issues we have stopped development of a spectral 
implementation of the level-set method. These problems arise because of properties 
of the function /, which is merely a tool to represent the actual surface of interest, 
/ = 0. This surface itself is well-behaved and smooth, suggesting it will be possible 
to evolve this surface directly. Geodesic and surface methods do precisely this, and so 
we focus on these two methods in the remainder of this paper. 

1 Even with re-initializations of /, as performed in [19] , the same argument applies to that time 
interval between re-initializations during which the topology of the EH changes. 
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2-4- Numerical Implementation 

Compared to the implementation of the geodesic method, implementing the surface 
method is somewhat more complex due to the presence of derivatives along the surface 
in (| 1 4 q|) . Apart from this, geodesic method and surface method share rather uniform 
implementation details. We shall first discuss those aspects that only apply to the 
surface method, and then follow with aspects applicable to both methods. 

We represent the surface r l (t,u,v) with spectral methods (e.g. 25J ) . These 
methods approximate a desired function t/(x, t) as a truncated expansion in basis 
functions (f>k, for instance Chebyshev polynomials or spherical harmonics: 

jV-l 

U(x,t)=J2Uk(t)<t>k(x), (17) 

fe=0 

where N is the order of the expansion. The fundamental advantage of spectral methods 
lies in their fast convergence: For smooth problems and a suitable choice of basis 
functions, the error of the approximation (I17|) decreases exponentially with the number 
of basis functions per dimension [25] . Derivatives of the function U are computed via 
the (analytically known) derivatives of the basis functions. Each set of basis functions 
has an associated set of collocation points x<; a matrix multiplication translates 
between function values at the collocation points, f/(xj), and spectral coefficients 
U k . 

For the surface method, we represent each Cartesian component of r l (t,u,v) (cf. 
([8])) as an expansion in scalar spherical harmonics, 

L rn=+l 

r%u,v) = E E ^mit) Y tm (u,v). (18) 

i=Q m=-l 

This expansion assumes that at fixed t, the surface has topology 52- Note that ([T5)l 
allows the surface to intersect itself, as necessary in a binary merger for t < icEH 
(cf. Figure [1]). Self- intersection is possible because the coordinates u and v are not 
assumed to be standard spherical angular coordinates, i.e. relations like cos(u) = 
zj \J x 2 + y 2 + z 2 will in general not hold. 

For spherical harmonics Y^ m (u, v), the collocation points form a rectangular grid 
in (u,v), with the u values chosen so that cos(u) are the roots of the Legendre 
polynomial of order L + 1, and with the v values being uniformly distributed in the 
interval [0,27r[. There are in total 

N = 2{L + 1) 2 (19) 

collocation points. The evolution equations (|13|) require derivatives d u r % and 
d v r % , which are computed by transformation to spectral coefficients, application of 
recurrence relations and inverse transform (using the SpherePack library [26)). These 
derivatives are then substituted into ([T3|) - (jl4c|) to compute ftr 1 , which is evolved at 
the collocation points. 

We represent each Cartesian component r % as an expansion in scalar spherical 
harmonics (see (|18p ) in order to re-use the infrastructure already developed for our 
spectral evolution code, which represents tensors of arbitrary rank in this manner to 
simplify our spectral expansions and to simplify communication of tensor quantities 
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across subdomains of different shapes (see, e.g., [371 HH])- An alternative approach 
would be to represent r l in terms of vector spherical harmonics, i.e., 

l m=+e 

r%u,v)=J2 E 4m(i)4(«.»)- ( 20 ) 
i=a m =-e 

The downside of choosing a scalar spherical harmonic representation is that the 
equation we impose on the highest order vector spherical harmonics is incorrect, and 
this leads to an instability. This difficulty with expanding vector quantities in a scalar 
spherical harmonic basis is cured |27| by performing the following "filtering" operation 
at each timestep: first transform r % to a vector spherical harmonic basis, then remove 
the I = L and £ = L — 1 coefficients, and then transform back. The removal of both 
the highest and second highest tensor harmonic modes is necessary, since transforming 
an n-th rank tensor from a tensor spherical harmonics to scalar spherical harmonics 
requires scalar harmonics of up to L scalar = l} ensor -f n . We filter two modes because 
we wish to correctly represent the spatial derivatives of r l (see (|14al) s ). which are 
effectively rank 2. 

The geodesic method simply evolves the ODEs (|7al) - ([76l) . While each geodesic 
is evolved independently, we find it nevertheless convenient to represent them as a 
two-dimensional grid, q i (t 1 u,v) where parameters u and v label each geodesic. We 
use the same parameters u and v for geodesic and surface method, and for this paper, 
we choose to locate the geodesies at the same (u, v) values as the collocation points 
of the surface method. We note that this choice is based on convenience to simplify 
comparison between the two methods; geodesies can be placed at any location, and 
indeed, we plan as a future upgrade of the geodesic method an adaptive placement of 
geodesies to help resolve interesting features like caustics. 

Let us now discuss aspects common to the implementation of the geodesic and 
surface methods: At some late time t — £ cn d long after merger, we initialize the EH 
surface by choosing it to be the AH at that time. Our AH finder parameterizes the 
radius of the AH as a function of standard azimuthal and longitudinal angles on 52, 
rAn(t e nd,0,<f>), i-e. 

(sin (9 cos \ 
sin 6 sin <p . (21) 
cos (9 / 

When initializing the event horizon surface, we choose (u, v) to coincide with the 
standard spherical angular coordinates (#,</>), i.e. we set 

r'fiend, u, v) = r AR (t cnd , u, v), surface method, (22) 
u, v) = r l AH (t cnd , u, v), geodesic method. (23) 

For the geodesic method we further set p l (t end , u, v) = s\hj where s^h is the unit 
normal to the apparent horizon, which is computed similarly to (|14a[) - (|14c[) . Time 
stepping is conducted using a 4th order Runge-Kutta algorithm. 

Both methods require interpolation of certain quantities like the spatial metric 7^ 
onto the grid points of the surface r*(i, u, v). For the spectral evolutions of the Caltech- 
Cornell group [20l EH [28] the evolution data is represented as spectral expansions in 
space (for each fixed time t) and spatial interpolation is performed by evaluating the 
appropriate spectral expansions (fTT[) at the desired spatial coordinates r % . Evolution 
data is available at discrete evolution times t n and temporal interpolation is performed 
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with 6-th order Lagrange interpolation (i.e. utilizing 3 time slices on either side ol the 
required time)0 

Finally, we define an area element \fh on the surface as the root of the determinant 
of the induced metric, 

h= _i_ detf^i^l ^Yi^A ■ ( 24 ) 

u v r o v ? J 

The area of the evolved surface is then given by 



A(t) = dA= / y/h{t,u,v) sin u du dv. (25) 



Explicitly pulling out the factor sin u in (|2"4"|) and (j2"5|) ensures that \fh is a constant 
for a coordinate sphere in Euclidean space; this will simplify Figure [11] below. Since 
all the geodesies (or surface grid points) are on a Legendre-Gauss grid, we compute 
the derivatives in (|24|) spectrally, and we evaluate ([25]) by Legendre-Gauss quadrature. 
For binary black hole mergers before merger, we sometimes evaluate h based on finite- 
difference derivatives d u r l and d v r l . This is discussed in detail in Section l4~4l 

3. Application to Kerr spacetime 

Initial tests of the event horizon finder are conducted using the Kerr spacetime in 
Kerr-Schild coordinates (See §33.6 of [29]): 

9^ = Vnu + ZHl^lv. (26) 

Here H is a scalar function of the coordinates, i]^ is the Minkowski metric, and is 
a null vector. In Cartesian coordinates (t,x,y,z), the functions H and / M for a black 
hole of mass M and dimensionless spin parameter a/M in the z direction are 

H - ~i — ; — 2~~2 > ( 27a > 

, _ A xr B L + ay yr B h - ax z \ 

where rBh(x,y, z) is the Boyer-Lindquist radial coordinate, defined by 

r| L = ^(x 2 +y 2 + z 2 - a 2 ) + (l (x 2 + y 2 + z 2 a 2 f + a 2 z 2 ^ ^ . (28) 

If we define the Kerr-Schild spherical coordinates in the straightforward way (r = 
x 2 + y 2 + z 2 , cos(8) = z/r, etc), we find that the event horizon of the Kerr black 
hole in these coordinates is given by 



/ r ~\~ r 2 a 2 

r K err(M) = W 2 + 2 + 2fl , (29) 

V r± + a z cos z 



where r + = M + y/M 2 — a 2 . The surface area of the event horizon is given by 



A KelI = 8ttM(M + y/M 2 -a 2 ). (30) 

+ The spectral spatial interpolation is computationally more expensive than temporal Lagrangian 
interpolation. Whenever the domain decomposition for the Einstein evolution is identical for all 
timesteps involved in a temporal interpolation, the time interpolation is performed before the 
spatial interpolation. In that case, only one spectral spatial interpolation is necessary (on the time- 
interpolated data), rather than six, speeding up the computation. 
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For our tests on Kerr spacetime we choose the same initial surface for both the 
surface and geodesic methods: a coordinate sphere of radius r = 2.5M, which does not 
coincide with the horizon. The evolution begins at t en< j = and proceeds backward 
in time towards negative t. Because we choose to place geodesies coincident with 
the collocation points of the surface method (see Section I2.4[) , we can use the highest 
angular index L as a measure of resolution. The total number of geodesies or grid 
points is given by (|19p . The choice of spin in the z direction is for convenience. We 
have repeated the numerical tests below for spins of several different orientations, and 
we find no substantial difference in cither stability or accuracy. 

In order to test our methods of finding an EH, we use two measures of error. The 
first measures the error in the coordinate location of the event horizon. We define 

Ar(u, v) = r(u, v) - r Kerr (0(u, v),<j)(u, v)). (31) 

where r(u,v), 6(u,v), and 4>(u,v) are the Kerr-Schild radial and angular 
coordinates of the surface, which are found from either the surface-method variables 
r l (u,v) = [x(u, v),y(u, v),z(u, v)] or the geodesic-method variables q l (u,v) = 
[x(u, v), y(u, v), z(u, v)] in the usual way, e.g., x(u,v) = r(u,v)sm9(u,v)cos(f)(u,v). 
Specifically, we will use the root-mean-square of Ar over all grid points or geodesies, 
which we shall denote by ||Ar||, as a global measure of the error. 

Our second error measure is the deviation of the area of our surface from the Kerr 
value, 

AA = A{t) - A K err, (32) 

where A(t) is determined by Equation (|25|) . 

Figure [3] shows errors in the AH surface as computed using the geodesic method 
for a Kerr black hole. The error measure ||Ar||, (f3Tj) . does not change with L because 
the evolution of each geodesic is independent of the total number of geodesies. The 
error measure \AA\, (|32|) . does depend on L, but only because the computation of 
the surface area depends on all geodesies. It is clear from Figure [3] that the geodesic 
method can stably model Kerr black holes of any spin. 

At i nd = 0, we start the EH finder with an initial surface that does not coincide 
with the EH of Kerr. Therefore, Figure [3] shows initial transients as the surface being 
followed by the EH finder approaches the EH of Kerr. Figure 2] shows an enlargement 
of this phase. We find that the tracked surface approaches the Kerr EH exponentially 
when integrating backward in time, 

||Ar|| oce* /T . (33) 

The time scale r depends on the spin of the Kerr background. It has been shown in 
a number of coordinate systems [TH US1 [S] that the e- folding time for a non-spinning 
black hole is r = AM. This is not true in all coordinate systems: for example, in 
Schwarzschild coordinates r = 2M. In | Appendix B[ we generalize this result to show 
that null geodesies, perturbed from the Kerr EH, diverge from the EH exponentially 
with an e- folding time equal to l/g#, where 



VM 2 - a 2 

9h = r 1 ; (34) 

2M (M + VM 2 - a 2 ) V ' 

is the surface gravity of the horizon in Kerr-Schild coordinates. In TableQ] we compare 
the numerically computed e- folding time r (obtained by least-squares fits) to gn, and 
find excellent agreement. 



Revisiting Event Horizon Finders 



13 



10 

B io" 3 

5 io- 6 



— 10 



10 



12 



io l 



10" 



< , 

— 10 



10 



10 



12 



15 



a/M = 0.6 
L = 5 



11 



a/M = 0.6 



L = 5,7,9,11,13 



3000 -2000 -1000 

t/M 



L= 11 

a/M = 0.99 







0.9 










0.8 
0.6 






0.4 N| 







i i , oi 



L= 11 



a/M = 0.99 



0.4 0.6 0.8 0.9 

v \ ^ 



0.2 



-2000 -1000 
t/M 



Figure 3. Geodesic method applied to a Kerr black hole. The top panels 

show the area difference between the computed and exact solution, normalized by 
the area of the exact solution. The bottom panels show the difference between 
the computed and exact location of the EH, as measured by (|31fl . These data are 
shown for two series of runs: In the left panels we keep the dimensionless spin 
of the black hole fixed at a/M = 0.6 and vary the resolution L of the EH finder. 
In the right panels we vary the spin parameter a/M at fixed resolution. In all 
cases, the EH finder starts at t = and the geodesies are evolved backward in 
time. 



Table 1. Exponential approach of the null surface to the correct event horizon 
location. Mgu represents the (dimensionless) surface-gravity of a Kerr black hole 
with spin a/M. M/t is the numerical rate of approach as determined by fits to 
the data shown in Figure [4] 
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Figure 4. Approach of the tracked null surface onto the event horizon of Kerr 
black holes with various spins. The symbols show the numerical data (the same 
data as in the lower right panel of Figure [3)1, and the solid lines are representative 
least-squares fits. Table [T] compares the numerically computed e-folding time to 
the surface gravity of the black hole. 



We now turn our attention to the surface method. For a Schwarzschild black hole, 
the surface method with the standard tensor spherical harmonic filtering is stable, as 
shown by the "F=0" line in the left panel of FigureO However, the method is unstable 
for spinning black holes and fails within about 10M for spin a/M = 0.6 (see the "F=0" 
line in the right panel of Figure O . 

Therefore, we perform additional filtering for spinning black holes. After each 
timestep, we compute 

R(u,v) = \J 5ijf l {u, v)ri (u, v), (35) 

expand R(u, v) in scalar spherical harmonics, 

l i 

7 / I 

(u,v), (36) 

1=0 m=-l 

and truncate the highest F modes of this expansion: 

R tm -► 0, for £> L-F. (37) 

From these filtered coefficients, we reconstruct the filtered radius-function Rp(u,v) 
and replace 

r 4 - (38) 

The right panel shows that with appropriate choice of F, the horizon of a Kerr 
black hole with spin a/M = 0.6 can be followed for thousands of M. Unfortunately, 
we do not understand the effect of F on stability, and therefore a parameter search 
through possible values for F is required. 

With this additional filtering in place, we now examine the convergence and 
accuracy of the surface method. Figure [5] shows the convergence behaviour of the 
surface method. From the top plots, we can see that for a black hole of moderate 
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t/M t/M 

Figure 5. Effect of filtering using J37H and (13S I t for the surface method. Shown 
are evolutions with the same angular resolution L = 18, but for different numbers 
F of truncated modes in 137> . Left panel: For a Schwarzschild black hole, the 
surface method is stable with or without this additional filtering. Right panel: 
For a Kerr black hole with a/M = 0.6, F = 7 performs best. The EH finder starts 
at t = and the surface is evolved backward in time. 



spin (a/M = 0.6), the surface method is accurate and convergent, although long-term 
stability issues remain. Also, the surface area computed by the surface method appears 
to be more accurate than the location of the surface, cf. upper vs. lower panels of 
Figure[Hl This arises, because for a small change SA l lm in an expansion coefficient A\ m 
in (fl"8"|) with I 7^ 0, the change in ||Ar|| is linear in 5 A) , whereas the change in area 
is quadratic. The high accuracy of Aeh is a welcome feature, since the EH area is one 
of the most important results of an EH finder. Unfortunately, the surface method is 
not capable of tracking the horizon for spins a/M > 0.8 for a useful length of time. 

While the geodesic method appears superior in these Kerr tests, there are two 
main benefits to implementing the surface method. Firstly, it is computationally more 
efficient. The bulk of processing time is spent on interpolating the metric data from the 
simulation, and the surface method requires the metric only (10 components) whereas 
the geodesic method requires the metric, as well as its spatial and time derivatives 
(50 components). Secondly, the surface method can be used to check the errors in the 
geodesic method in circumstances where the surface method performs well, i.e. lower 
spins. 

For these tests, the initial set of geodesies (or surface) is chosen to be a sphere of 
radius 2.5M. In this case it requires a time > 100M for either method to converge onto 
the actual event horizon. This shows that for cases in which the actual EH is unknown, 
it is important to have a near-stationary situation at the end of the simulation, so that 
the initial guess (generally taken to be the AH) has time to converge onto the true 
EH. The length of this interval will depend on the desired accuracy, the quality of the 
initial guess and the spin of the black hole. For example, during a time At = 10/ gn 
(i.e. 40M for a/M = 0, but 160Af for a/M = 0.99) the tracked surface will have 
approached the EH to a fraction e~ 10 ~ 5 • 10~ 5 of the distance between the initial 
guess and the EH. 
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Figure 6. Surface Method, applied to a Kerr black hole. The top panels show 
the normalized area difference between the computed and exact solution. The 
bottom panels show the difference between the computed and exact location of 
the EH, as measured by H31I I. These data are shown for two series of runs: In the 
left panels we keep the dimensionless spin of the black hole fixed at a/M = 0.6 
and vary the resolution L of the EH finder. In the right panels we vary the spin 
a at fixed resolution. The value F denotes the number of truncated modes during 
filtering according to H37I I. For each case, we show the value of F that provides 
the most accurate evolution. Also, in all cases, the EH finder starts at i cnt j = 
and the surface is evolved backward in time. Compare to Figure [3] 



4. Head-on Binary Black Hole Merger 

4-1. Details of BBH evolution 

When looking for a straightforward dynamical spacetime where tracking the event 
horizon is of interest, one of the standard scenarios is the head on merger of two 
equal-mass non-spinning black holes [TTJ [TH [TJl E3 HSj- First, the SpEC code is 
utilized to evolve the solution of Einstein's equations for the head-on merger. Initially 
the holes are at rest, r ~ 4.5M apart, where M = Ma + Mb is the total mass at t = 
(because the black holes are non-spinning, we take the irreducible mass as the black 
hole mass, M A / B = M in . A / B = ^Maha/b/(167t)). Initial data is constructed by 
solving the conformal thin sandwich equations [301 EI] with the same setup as in [28] , 
but setting the orbital frequency SIq = 0. This data is then evolved with the SpEC 
code using the dual coordinate frame technique described in [55] and with a domain 
decomposition with two excision spheres. A common apparent horizon appears at 
t = tcAH = 17.83M. Shortly thereafter, at t — t re grid = 18.96M, the original domain 
decomposition with two excision boundaries is replaced by a set of concentric spherical 
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Figure 7. Evolution of a head-on BBH merger: normalized constraint 

violations. The left panel shows the complete evolution. The right panel enlarges 
the time around merger, with formation of a common apparent horizon and time 
of regridding indicated by 'CAH' and 'regrid', respectively. The discontinuity at 
iCAH arises because the constraints are computed only outside the common AH 



for t > *cah- At t. 



regrid 



the constraints jump because of the different numerical 



truncation error of the ringdown domain decomposition. 



shells with one larger excision boundary. The new excision boundary lies somewhat 
inside the common apparent horizon, but outside the original excision boundaries. The 
region very close to the original excision boundaries, and between them, is dropped, 
and is no longer evolved. Data is interpolated from the highest resolution merger run 
onto three resolutions of this new domain decomposition. The simulation is continued 
up to t — 95M and the final mass of the merged black hole is Mfi na i = 0.9493M. 

The simulation is performed at three progressively higher resolutions, named 'NO' 
through 'N2'. The SpEC code does not strictly enforce the Hamiltonian or momentum 
constraints, nor the artificial constraints that arise from the first-order reduction of 
the Generalized Harmonic formulation of Einstein's equations [32] ■ As such, it is 
important to monitor the values of these constraints during the simulation, as shown 
in Fig. [71 We normalize the constraints by an appropriate norm of the derivatives of 
the evolved variables (see (71) of [35] for the precise definition) and integrate constraint 
violations and normalization only outside the two individual apparent horizons or the 
common apparent horizon for this run. 

4-2. EH finder behaviour 

Since the EH finder follows the EH backward in time, we begin our discussion with 
the ringdown phase of the head-on merger. Initial data for both the geodesic and 
surface methods is taken from the apparent horizon at t = 81.24M, about 60M after 
appearance of a common AH. 

We run both the geodesic and surface methods for angular resolutions L = 
7, 15, 23, . . . , 47 and compute the area A(t) of the tracked surface for these runs. We 
do not employ filtering as per ([3"7| for the surface method. 

Figure [8] plots the relative differences between A(t) computed with different 
angular resolution. This plot exhibits several noteworthy features, which we discuss 
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Figure 8. Effect of changing the resolution of the EH finder when applied to 
the BBH evolution at fixed high resolution. Shown are relative differences in the 
area A(t) of the tracked surface. The label "7-15" denotes the difference between 
simulations with L\ = 7 and L2 = 15, normalized by A(t) of Li- Vertical lines 
on the graph denote the formation of common event and apparent horizons. Note 
that the time scale of both plots change at t/M = 25. 



in the next few paragraphs: 

During the ringdown phase, t > 20M, both the surface and geodesic methods 
perform admirably: Even at low resolution L — 7, the area is computed to better 
than 10 -6 and this error drops rapidly below 10~ 12 as L is increased. The rapid 
convergence with L in the ringdown regime is not too surprising, because the angular 
resolution of the merger simulation is Evolution = 25. Therefore, angular modes t > 25 
of the EH finder carry only information about the way in which the surface parameters 
(u, v) deviate from the 4>) coordinates of the simulation. As can be seen from the 
excellent convergence for t > 20M in Figure[8j such deviations are not very important. 
We also note that the long-term instability exhibited by the surface method during 
the Kerr test is not apparent. 

Close to merger and before merger, t < 20M, the tracked surface becomes very 
distorted and therefore requires much higher angular resolution. This is apparent in 
the comparatively larger errors in A{t) for icEH < t < 20M . In this time interval, the 
errors in the surface method grow more rapidly than those of the geodesic method. We 
attribute this to a degradation of the convergence rate of the spectral expansion [)18p . 
The surface method relies on the spectral expansion in an essential way to compute 
the derivatives that enter into (|14aj) . In contrast, evolution of geodesies is independent 
of the spectral expansion and the spectral series is used only to compute the surface 
area via (1251) . 

At the point of merger, when the surface being tracked by the EH finders intersects 
itself for the first time, the error in the area-computation suddenly increases drastically 
in either method. The reasons for this are quite different for the two methods: The 
geodesic method evolves individual geodesies perfectly fine through icEH- The large 
errors in Figure [8] arise because of the use of spectral integration to compute the surface 
area: At a caustic, the surface-area element y/h, tends to zero, resulting in a 
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Figure 9. Error estimates for the surface and geodesic methods, with surface 
resolution L = 47. The left panels show the root-mean-square pointwise deviation 
between the different runs, whereas the right panels show the differences in the 
surface area. The lines labelled "G:N#-N#" ("S:N#:-N#") in the upper panels 
show the difference between the geodesic method (surface method) when applied 
to merger simulations of different resolution N. The lines labelled "N#:S-G" in 
the lower panels show the differences between the surface and geodesic methods 
for a given N (where NO, Nl, and N2 are resolutions of the merger simulation). 
Note that the time scale of all plots change at t = 25M. 

non-smooth integrand in the area integral (|25[) , destroying exponential convergence of 
the spectral area integration. Below, we will explain how we employ finite-difference 
integration instead. We shall address area calculation for t < £ceh in Section 14.41 
where we also discuss how to compute the area of the EH excluding the future 
generators of the EH. 

The surface method exhibits additional, more fundamental, problems at icEH, 
when the surface being tracked intersects itself in a caustic with \fh — > 0. At such 
a point, the tangents to the surface, d u r l and d v r l are either no longer linearly 
independent, or one of them is zero, cf. ([24]) . Therefore the surface normal s l in 
(fT4"q[) is ill-defined. 
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While the surface method presently cannot evolve through merger, it nevertheless 
yields valuable consistency checks with the geodesic method during the ringdown 
phase. Figure [5] presents such a comparison between the two methods and examines 
the effect of varying the resolution of the underlying binary black hole simulation. 
The top panels show differences between the results of the geodesic method applied 
to evolutions with different resolutions (labelled "G:N#-N#"). As the underlying 
resolution is increased, the differences become smaller. Likewise, the lines labelled 
"S:N#-N#" show the analogous differences when running the surface method. When 
the surface method works, t > 15M, it is more accurate than the geodesic method. 
For times close to the formation of the common event horizon, t < f 5M, errors in the 
surface method grow very rapidly and render our current implementation essentially 
useless. The bottom panels of Fig. [9] show differences between surface and geodesic 
method at the same resolution of the evolved data. This difference decreases with 
increasing N, as it should. During ringdown, t > 15M, the difference is essentially 
equal to the error in the geodesic method; for t < 15 M it is dominated by errors in 
the surface method. 

The right panels in Figure [9] examine the surface area A(t). No clear convergence 
is apparent for t > 20M, perhaps because the surface area of the event horizon can 
be calculated with great accuracy even at low values of N. Given the lack of clear 
convergence, we shall take as our error estimate for the post-merger area the square 
sum of the following three error measures: a) the change in A(t) between the geodesic 
method applied to the head-on simulation at the two highest resolutions (i.e. "G:N1- 
N2"), b) the change in A(t) between the geodesic and surface methods (i.e. "N2:S-G") 
and finally c) the change in A(t) in the geodesic method at L = 47, N2 when doubling 
the timestep (from 0.056M to Q.112M; the effect of this is small and not shown in 
Figure [9]). This combined error estimate is plotted in Figure [TOl 

4-3. Quasinormal Modes during Ringdown 

After the merger, the distorted merged black hole rings down into a stationary black 
hole. During this phase, the area of the event horizon, Aeh will approach its final 
value Api na i, and one expects that the apparent horizon approaches the event horizon. 
This is explored in Figure [TDJ This plot also contains the error estimates obtained 
from Figures [8] and [9l Figure [TO] shows that the areas of the common AH and EH 
differ by about f 0% when the common AH first appears, though this difference drops 
to 0.1% within about 3M. After this rapid initial drop, the ringdown is clearly 
apparent. The area of both EH and AH approach their final area exponentially, 
and this approach is resolved through about five orders of magnitude. A least- 
squares fit of log [Af — Aeh(£)] to the function C — X h s t for 30M < t < 70M, yields 
A bs = There are furthermore periodic features visible in the EH and AH 

areas, with seven periods clearly distinguishable. The period of oscillation is found to 
be t osc = 8.00A/, therefore uj ohs = 0.745Mg n 1 al . 

Decay rate A b s and frequency w b s can be related to quasi-normal modes of a 
Schwarzschild black hole as follows: The quasinormal mode parameters of a perturbed 
black hole are typically defined with reference to oscillations in the metric fields, which 
can be written as 

Sg^ oc e~ xt sin(u;t), (39) 

where A is the decay coefficient and u) is the angular frequency of the metric oscillation. 




Figure 10. Surface area differences between the EH, AH, and the final area, 
normalized by the final area. Also plotted are error estimates for A^n(t) and 
A AH (t). 



Therefore 

Sij^n oc — \e~ xt sin{u)t) + uje~ xt cos(Lut). (40) 

The energy flux through the horizon, and therefore the change of its mass is M oc 
|<%«/| 2 ! so we have 

A e~ 2Xt 

— oc M oc [A 2 + uj 2 + (lu 2 - \ 2 )cos(2iot) - Xcjsin(2ut)]. (41) 

J\ 2 

Thus the observed values (A b s , w b s ) should be twice the values (A,w) of a quasi- 
normal mode. Indeed, the lowest quasinormal mode of a perturbed Schwarzschild 
black hole is the t = 2, n = mode, with [33J A 20 = 0.08896^/,^ and uj 2 q = 
0.37367]^^. Consistent with flU]), we find that A obs - 2A 20 = O.OOSAffT 1 ^, and 
co obs - 2uj 20 = 0.002^4"^. 

4- 4- Treatment of Merger 

Before examining the merger phase in detail, we must develop tools to analyse the 
topology change the event horizon undergoes during merger. As seen in Figure [TJ 
prior to merger, the surface found by the event horizon finder is the union of the two 
individual event horizons and the set of future generators of the joint event horizon. 
The event horizon itself consists of two topological spheres. At merger, t — £ceh, 
the topology of the event horizon changes to a sphere. For t < tcEH, generators 
of the event horizon continuously enter the event horizon at the cusps on the event 
horizons of the two approaching holes. The geodesic method traces geodesies perfectly 
fine through merger back to the start of the head-on binary black-hole evolution, and 
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Figure 11. Left panel: Area element yh along a few representative geodesies 
during the head-on merger. Each geodesic is labelled by the angle a between the 
initial location of the geodesic (at t en d) and the axis of symmetry. Three types 
of behaviour are apparent: Geodesies entering the horizon from Z~ (a = 85°); 
geodesies entering the horizon from an area in the vicinity of the individual event 
horizons before merger (a = 72° or 76°), and geodesies remaining on the horizon 
throughout. The right panels show the time derivative of yh, highlighting the 
clear signature when a geodesic enters the horizon. 



the trajectories of the geodesies are convergent as the resolution of the underlying 
evolution is increased, see the left panels of Figure O In this section, we address two 
questions relevant to analysing the output of the geodesic method: First, when going 
toward earlier times, some geodesies leave the event horizon; how does one decide 
whether a given geodesic is still on the event horizon, or whether it is merely a future 
generator of the event horizon? Second, how can one compute the area of the event 
horizon (i.e. not counting the area of the locus of future generators)? 

Let us first consider the area element yfh of the EH surface, with h given by (|24jl . 
which requires derivatives d u , d v along the surface, thus connecting neighbouring 
geodesies. Because we place the geodesies at a (u,v) grid consistent with spherical 
harmonic basis functions, we can use spectral differentiation to compute these 
derivatives (and have done so, up to this point in the paper). Convergence of this 
spectral expansion, however, becomes increasingly slow for t < £ceh, an d therefore, 
we compute henceforth the derivatives d u r l and d v r l with second order finite difference 
stencils. 

Figure [IT] plots the area element yh as a function of time for a few representative 
geodesies. This figure was obtained from our highest resolution run using 20,000 
geodesies. To reduce CPU cost, these geodesies were initialized at t = 19.8M from 
the L = 48 run of the surface method. For some geodesies in FigureQTJ \ h approaches 
zero at a certain time. This feature can be used to determine whether a given geodesic 
is still on the horizon: We first note that the change of area element along a given 
null geodesic (i.e. for fixed u, v) is proportional to the expansion of this particular 
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geodesic: 

9t lo g ^ = ^J) ocfl. (42) 

The constant of proportionality depends on the parameterization of the null geodesic. 
Note that by Raychaudhuri's equation, the expansion of a generator of the event 
horizon must be non- negative, 9 > 0. Figure QTJ shows the area element as a function 
of time for a few representative geodesies. 

At late time t = t nd where the final black hole has settled down, we start with 
geodesies on the apparent horizon, which will be very close to the event horizon. 
Therefore, we assume that at t on d all tracked geodesies are generators of the event 
horizon. Consistently with this assumption, Figure QT] shows that dt \ogV~h starts out 
very close to zero, and increases as we approach the dynamical time region around 
merger. If a generator remains on the event horizon, dt log y/h will eventually decrease 
again and approach zero at very early times before the merger. Generators leaving 
the event horizon must do so at points where generators cross, according to a theorem 
by Penrose 34. 29 . For the head-on merger, at such a point nearby geodesies cross 
and pass through each other. Just after the geodesic enters the horizon, the horizon 
generators diverge from each other and their expansion is positive (and so is d t log y/h) . 
Just before the caustic points, nearby future generators of the event horizon converge 
toward the caustic point with negative expansion. In fact, at the caustic, dty/h changes 
sign discontinuously, as can be seen in Figure [TTJ 

Therefore, the largest time at which the expansion of a geodesic passes through 
zero will be the time it joins the event horizon, 

>°; I < 43 > 

In practice, we keep track of (|43j) with a mask function /m(u,v), which is initially 
identical to unity. As we evolve backward in time, we evaluate dt log y/h at each time 
step, and if it drops below some tolerance — tol for a point (u , vq) we set /mIuq, vq) = 
for that geodesic. The tolerance tol is necessary to avoid misidentifications due to 
numerical truncation error at very early or late times, where dt log y/h — > for event 
horizon generators. Because dt\ogy/h changes so rapidly at a caustic, the precise 
value for tol is not very important; we use tol = 10 . 

For generic situations, generators can also leave the EH at points where 
finitely separated generators cross (a "cross-over point" in the language of Husa & 
Winicour [35]). At such points, y/h remains positive, and criterion (|43p reduces to a 
necessary but not sufficient condition that a generator has left the horizon, i.e. t- io \ n 
from (|43|) will be a lower bound for the actual time when a particular geodesic leaves 
the horizon. Cross-over points could be diagnosed by monitoring the minimal distance 
between every pair of followed geodesies, and we shall discuss this point in more detail 
in a future publication. 

The area of the event horizon (consisting of the two disjoint components for 
t < ^ceh) is found by multiplying y/h by the mask function /m, and integrating: 



j4eh = / /m {u, v) y/ h(u, v) sin u du dv. (44) 



For t < icEHj there are two major sources of error in this integral: First, each geodesic 
can either be on or off the horizon. When Jm changes discontinuously from 1 to for a 
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Figure 12. Convergence of the surface area of the event horizon during merger. 
The lower plot shows results for placement of the geodesic pole parallel to the axis 
of symmetry (i.e. consistent with axisymmetry) , the upper plot has a geodesic axis 
perpendicular to the axis of symmetry of the merger. In both cases, geodesies are 
tracked using the geodesic method; derivatives for \/h (cf. 12411 ) are computed 
with finite-differences; geodesies are removed from the event horizon based on 
(I4.'{1 . Lines are the difference between each resolution and the next highest. 



geodesic, the area of the event horizon will change discontinuously. Note that this will 
occur at different times for different resolutions. The severity of this effect will depend 
on how many geodesies enter the horizon simultaneously, as illustrated by Figure 1121 
This figure shows the convergence of the event horizon area with increasing number 
of geodesies, and for two distinct orientations of the geodesies. In either case, the 
geodesies are initialized at t = 19.86A/ from the L = 47 surface method determining 
the event horizon during ringdown, and in either case the geodesies are placed on a 
rectangular (u, v) grid as detailed in Section l2~4l In the lower panel of Figure [T2l the 
geodesies are oriented respecting the axisymmetry (i.e. the u = polar axis is aligned 
with the axis of symmetry), whereas in the upper panel the u = axis is perpendicular 
to the axis of symmetry. The lower panel of Figure [T^J with geodesies respecting the 
symmetry, shows much larger variations in the area as the resolution is increased. This 
arises because due to the symmetry, a full ring of geodesies leaves simultaneously, 
thus amplifying the discontinuity of A-En(t). For perpendicular orientation of the 
geodesies, individual geodesies leave the horizon, resulting in smaller jumps; this is 
the configuration we will use in the next section to examine the physics of the black 
hole merger. 

The second source of error in the evaluation of (f44|) arises because the integrand 
is not smooth once geodesies have left the horizon. For fixed t < £ceh , V~h approaches 
zero linearly toward the caustic; off the horizon, fuVh = by virtue of the 
mask function, so overall, the integrand is only continuous, and we cannot expect 
exponential convergence of the integral, despite using a Gauss-quadrature formula to 
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t/M 

Figure 13. Area of the event horizon and of the apparent horizons before merger 
and during merger. The vertical dotted lines indicate formation of a common 
event horizon and appearance of a common apparent horizon; the inset shows an 
enlargement for early time. 

evaluate 

4- 5. Analysis of Merger Phase 

When evolving geodesies backward, we find that the first geodesic leaves the horizon 
at tcEH — 14.58M, the time of merger. However, it should be noted that the point 
at which an observer sees the EH change topology is not invariant because the curve 
traced by the cusps of the two black holes is spacelike [36 . Figure [13] shows the 
surface area of the EH and the common and individual AHs during the merger phase. 
The common apparent horizon forms at icAH = 17.8M, and we track the individual 
apparent horizons up to t = 18. 8M . The area of the individual apparent horizons 
is remarkably constant. Up to formation of the common event horizon, its fractional 
increase is less than 10 , up to common apparent horizon, its fractional increase is 

5- 10 -5 , and even when we stop tracking the inner horizons, their area has increased by 
only 1.6 • 10~ 4 . In contrast, Aeh varies significantly more and at significantly earlier 
times, as can be seen from the inset. 

To examine the relation between individual apparent horizons and event horizons, 
we plot in Figure [14] the difference AA = A^h — (^ah,a + ^4ah,b)- For times 
6 < t/M < 14, AA grows exponentially with an e-folding time of 1.95M. This 
c-folding time is within a few percent of the surface gravity of a black hole with 
the initial mass of the black holes in the head-on simulation. This confirms that as 
geodesies are integrated backwards in time, the individual components of the event 
horizon approach the individual apparent horizons with the expected rate. If our code 

* For t > tcEHi ^EH m Figure [3~2l is limited by the finite-difference derivatives used to compute \fh. 
Better accuracy can be obtained using spectral derivatives, as can be seen from the right panels of 
Figure [9] For the analysis of the merger below, this difference is invisible. 
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were free from all numerical errors, the curve in Figure [H] would continue to decrease 
exponentially as one proceeds backwards in time. Instead, this curve saturates at 
AA/A as 0.1% at t = 0, and in addition, a feature in AA appears at t ~ 5M because 
the EH area falls below ^4ah,a + ^ah.b and therefore AA changes sign. These effects 
are due to numerical errors, particulary finite-difference errors in the computation of 
Vh and the use of a finite number of geodesies. If one wishes to achieve much better 
than 0.1% accuracy of the event horizon surface area at very early times when the 
two holes are widely separated, the EH must be split into two individual surfaces to 
be evolved separately. 



5. Conclusion 



This paper examines three different methods for locating event horizons in dynamical 
black hole spacetimes, the geodesic method, the surface method and the level-set 
method. All three methods rely on the principle that outgoing geodesies exponentially 
approach the event horizon when followed backward in time. We implement both the 
geodesic and surface methods, the latter one implemented without the assumption 
of axisymmetry as done in earlier work [12]. Overall, we find that the geodesic 
method is more robust, with the capability to accurately follow highly spinning black 
holes (tested up to a/M = 0.99), as well as the merger of two black holes. For the 
head-on merger, we find that the surface-area element \J~h of the geodesic congruence 
is an excellent diagnostic of whether and when a geodesic joins the event horizon, 
cf. (|43p . In more generic situations, this criterion might have to be amended by a 
second test for crossing of geodesies that are initially (at t — i C nd) separated by a 
finite amount. Errors due to tangential drift of the geodesies — as explained in [12"]— 
are not apparent in our simulations. The observed good properties of the geodesic 
method might be related to the improvements in accuracy of the spacetime metric 
since the early tests [T2], as well as the ability to interpolate the metric spectrally to 
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the geodesic locations. Because each geodesic is evolved independently, the geodesic 
method parallelizes trivially. Tracking of the cusp of the disjoint components of the 
event horizon before merger, as well as computation of Aeh is currently not highly 
accurate, as comparatively few geodesies cover the region close to the cusps. Our 
current scheme switches from the surface method to a large number of geodesies some 
time after merger where the surface method is still very accurate, say at to- We plan 
on adaptively placing additional geodesies at to based on where cusps occur. 

The surface method is less robust and exhibits a long-term instability when 
applied to Kerr black holes with spins a/M < 0.6, and rapid blow-up for larger spins. 
Nevertheless during the ringdown phase t > icEH of the head-on merger, the surface 
method locates the event horizon with comparable accuracy to the geodesic method 
and provides an important independent test of the geodesic method. However, when 
the surface being tracked self-intersects in a caustic point, our current method for 
defining the normal breaks down because dr % jdv = in p4al) - (|14cj) . and thus our 
current implementation of the surface method fails. 

The level-set method, finally, is not implemented in this paper. It requires 
boundary conditions for the level-set function /; furthermore / can become singular 
during a black hole merger. Both reasons made it unduly difficult to implement this 
method in our spectral code. In conclusion, we find that the geodesic method, the 
oldest of the three methods considered, to be the most accurate and useful in our 
tests. 

Turning our attention to applications of the event horizon finders, Figure |4] 
presents a new quantitative test of event horizon finders: When finding the EH of 
a Kerr black hole starting away from the true horizon, does the tracked null surface 
approach the true event horizon with the correct rate, namely the surface gravity 
gn7 Table Q] confirms this for the geodesic method. For the head-on merger, both 
geodesic and surface method perform admirably during the ringdown phase, where we 
are able to clearly observe the quasinormal ringing of the single merged black hole. 
For both the event and apparent horizons, the frequency and damping time of the 
ringing matches the (£ = 2,n = 0) mode of the Schwarzschild quasinormal ringing 
spectrum to within 2% for the decay rate and 0.3% for the frequency. 

Furthermore, we find that the apparent horizons provide an excellent 
approximation to the event horizon for the head-on merger very early before the 
merger, and very late after the ringdown. Thus, while in principle the apparent 
horizon is slice-dependent and there is no guarantee that it should coincide with the 
event horizon, in practice no such behaviour is found. 

Finally, perhaps surprisingly, for the head-on binary black hole merger only some 
of the future null generators of the horizon start at past null infinity. A significant 
fraction of the generators rather start close to the individual event horizons of the 
black holes before merger. This can be seen in the spacetime diagram in Figure [15j 
most clearly for the geodesic pointed to with an arrow. These geodesies begin to 
diverge from the individual event horizon as the second black hole approaches. The 
increased gravity of both black holes causes such geodesies then to "turn around" and 
join the event horizon at the seam of the pair of pants. 

In the future we plan to study event horizons in the more diverse black hole 
scenarios currently being simulated: mergers of inspiraling black holes; spinning 
and/or non-equal mass binary black holes, as well as black hole Neutron star mergers. 
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Appendix A. Points in the Surface Method follow geodesies 

Consider a 2-dimensional family of null geodesies, q^{t, it, v), where u, v label different 
geodesies. Assume the parameter t along the geodesic coincides with the coordinate 
time of the underlying black hole simulation, i.e. q°(t,u,v) = t. This family of 
geodesies traces out a three-dimensional null surface Af, parameterized by coordinates 
t,u,v: it, u), where t is the parameter along each null curve, and it, v are the 

parameters relating each null curve to nearby null curves. In this parameterization, 
we can write the outgoing null normal = dq t */dt\ u , i.e. a coordinate derivative 

£ = dt within the (t, u, v) coordinates of Af. Displacement vectors that relate each 
null curve to its neighbours are given by rh = d/du, n = d/dv. Since coordinate 
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derivatives commute, we have 

fV p ro" = mrVpF. (A.l) 

Let us consider the rate of change of the inner product £ M m M as we change the time t 
along a geodesic (i.e. for fixed u and v): 

Wrap) = fV„(fm M ) = m/V„(f ) + f/V„(m"). (A.2) 
From IjA.ip . the second term of (|A.2|) vanishes, 

^VvK) = ^m"V„(f) = im"V„(*%) = 0. (A.3) 

Substituting the formula for parallel transport of along the geodesies, ' v £^ = k£^ 
(with k = if t is affine), (IA.2P finally becomes 

fit(*"nv) = m/V^f) = «m^". (A.4) 
A similar calculation results in dt{(-^n v ) = nl^n^. 

So far, this appendix only discusses the geodesic method. We now use the results 
just obtained to show that surface and geodesic methods will construct the same null 
surface J\f. Both methods start with the same two-dimensional surface at some late 
time to, and the tangent <p(ioi u, v) to the geodesies at to is chosen to be normal to the 
2-surface. Therefore, at to, ^ =<?'*, and the surfaces resulting from evolving both the 
geodesic and surface methods will coincide at times infinitesimally near to- Because 
i^rrifi = t^n^ = initially, (1A.4|) implies that € M m M = l^n^ — at all other times. 
Thus, the tangent to the geodesies always remains orthogonal to the surface described 
by the positions of all the geodesies at a given time t. Since <p is normal to that 
surface, null, outgoing, and has q° = 1, it is identical at all times to as constructed 
by the surface method. Therefore, we see that the surfaces obtained by the geodesic 
and surface methods agree, and both techniques trace out the same Af given the same 
initial conditions. 



Appendix B. Proof of Surface Gravity conjecture 

We consider a null geodesic <z M (i) that asymptotes to a horizon generator qg(t) for 
t — > — oo, i.e. 

9*(t) = + (B.l) 

with Sq^(t) — * as t — > -co. In the discussion of Figured] we have asserted that 

Sq^(t) cx e 9Ht , (B.2) 

where gn is the surface gravity of the black hole, and where the coordinates are 
Kerr-Schild coordinates, cf. [|26 p -(|28 p . To confirm this assertion, one can substitute 
(jB.ip into the geodesic equation and expand to linear order in 8q^ (where we assume 
that 5q^ , 5q^ , and 5q^ are of the same order). One then needs to show that the 
resulting linear equation indeed has the solution (|B.2p . 

The linearization of the geodesic equation is most easily performed in adopted 
coordinates. We have performed the analysis in "rotating spheroidal Kerr-Schild 
coordinates" x M = (t, Tbl , 0, <t>) , related to the standard Kerr-Schild coordinates of 



(|26l) - (|2"8p by the coordinate transformation 

x = J r| L + a 2 sin 9 cos (4> + Vint) , (B.3) 

V = \Ai L + fl2 sin 6 sin O + n H t) , (B.4) 

z = tbl cos0. (B.5) 
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The time t is not transformed. Horizon generators have the form = [t, r + , 6q, 4>q\, 
with r + = M + V M 2 — a 2 and 9o,4>q constants, i.e. q^ oc [1,0,0,0]. In these 
coordinates, we have considered the geodesic equation in affine parameterization, ([2]) 
and have indeed confirmed 

Sq"' cx e 9Ht (B.6) 

to leading order in 8q^ . Exponential divergence from a horizon generator — as in 
(|B.6p — is a property present in a quite general class of coordinate systems. For 
instance, consider the coordinate transformation 

t' = t + f(x l ), x 4 ' = x l '(x l ), (B.7) 

where the Jacobian dx % jdx % and its inverse are finite in a neighborhood of the horizon. 
In this case, 8q^ and 8q^ are related merely by a multiplication by the Jacobian, so 
the exponential behavior e 9nt is the same in both coordinate systems. The coordinate 
transformation (|B.3[) - (|B.5[) falls into this class, and therefore (IB. 61) implies (|B.2|) . 
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